#include <iostream>
#include <cmath>
#include "/home/mjarvis/misc/nr/ran0.cpp"

const int NITER = 3;
const int MAXCOUNT = 10000;

const double omm_min[NITER] = {0.15, 0.05, 0.0};
const double omm_max[NITER] = {0.4, 0.5, 1.0};
const double omb_min[NITER] = {0.01, 0.0, 0.0};
const double omb_max[NITER] = {0.04, 0.05, 0.2};
const double w_min[NITER] = {-1.2, -2.0, -3.0};
const double w_max[NITER] = {-0.6, -0.3, 0.0};
const double A_min[NITER] = {0.4, 0.1, 0.0};
const double A_max[NITER] = {1.1, 2.0, 4.0};
const double n_min[NITER] = {0.9, 0.8, 0.5};
const double n_max[NITER] = {1.05, 1.2, 1.5};
const double h_min[NITER] = {0.6, 0.5, 0.3};
const double h_max[NITER] = {0.8, 0.9, 1.0};

double range_rand(double min, double max)
{
  return min + (max-min)*ran0();
}

int main()
{
  double vol0 = 0.;

  for(int iter = 0; iter < NITER; iter++) {
    double vol = (omm_max[iter]-omm_min[iter]) * (h_max[iter]-h_min[iter]) *
      (omb_max[iter]-omb_min[iter]) * (w_max[iter]-w_min[iter]) *
      (A_max[iter]-A_min[iter]) * (n_max[iter]-n_min[iter]);
    if (iter == 0) { vol0 = vol; }
    double volratio = vol/vol0;

    for(int count = 0; count < MAXCOUNT; count++) {
      double omm = range_rand(omm_min[iter],omm_max[iter]);
      double h = range_rand(h_min[iter],h_max[iter]);

      double Omega = omm/h/h;
      if (Omega > 1.) { count--; continue; }

      double omb = range_rand(omb_min[iter],omb_max[iter]);
      double w = range_rand(w_min[iter],w_max[iter]);
      double A = range_rand(A_min[iter],A_max[iter]);
      double n = range_rand(n_min[iter],n_max[iter]);

      std::cout << omb <<"  "<< omm <<"  0.  "<< n <<"  0.  ";
      std::cout << w <<"  "<< A <<"  "<<h <<"  "<<volratio<<std::endl;
    }
  }
}
